Conceptual Model
Contents
# import some utilities
# definition of some plotting tools
%matplotlib inline
from bokeh.layouts import gridplot
from bokeh.layouts import layout
from bokeh.plotting import figure, output_notebook, show
from bokeh.models import ColumnDataSource
from bokeh.plotting import output_file
from bokeh.palettes import Category10
from bokeh.io import export_png
from bokeh.models import Legend
from bokeh.models import Range1d
from bokeh.models import ColorBar
from bokeh.models import LinearAxis
from bokeh.models import DatetimeTickFormatter
# select a palette
#from bokeh.palettes import Dark2_25 as palette
from bokeh.palettes import Spectral11, Category20b
from bokeh.palettes import Spectral6
from bokeh.plotting import output_file
from bokeh.transform import linear_cmap
import itertools
import datetime
import itertools
import numpy
import os
import re
import sys
from netCDF4 import Dataset
import pandas
import logging
import math
def color_gen():
yield from itertools.cycle(Category10[10])
color = color_gen()
from bokeh.resources import INLINE
output_notebook(INLINE)
export = False
import logging
import sys
sys.path.append(os.path.join(os.environ["HOME"],"Projects","SeaFlect","python"))
from DataPoolAccessModule import AeolusDataClass
from DataPoolAccessModule import QualityControlClass
from SurfaceReflectanceModule import SurfaceReflectanceModelClass
def make_plot(title, hist, edges):
p = figure(title=title, tools='', background_fill_color="#fafafa")
p.quad(top=hist, bottom=0, left=edges[:-1], right=edges[1:],
fill_color="navy", line_color="white", alpha=0.5)
p.y_range.start = 0
#p.legend.location = "center_right"
#p.legend.background_fill_color = "#fefefe"
p.xaxis.axis_label = 'x'
p.yaxis.axis_label = 'Pr(x)'
p.grid.grid_line_color="white"
return p
verbose_level = 'info'
logger=logging.getLogger()
logging.basicConfig(level=getattr(logging, verbose_level.upper()))
Cm = {}
Li = {}
wind_speed = numpy.linspace(2, 24, 23)
theory=SurfaceReflectanceModelClass(logger=logger)
theory.Ro_fixed=0.016
Li["low"] = numpy.array([theory.surface_reflectance(model="li",u=key, theta=33) for key in wind_speed]).transpose()
Cm["low"] = numpy.array([theory.surface_reflectance(model="concept",u=key, theta=33) for key in wind_speed]).transpose()
theory.Ro_fixed=0.08
Li["high"] = numpy.array([theory.surface_reflectance(model="li",u=key, theta=33) for key in wind_speed]).transpose()
Cm["high"] = numpy.array([theory.surface_reflectance(model="concept",u=key, theta=33) for key in wind_speed]).transpose()
AeolusProductObject = AeolusDataClass(debug=False, \
generate=False,
logger=logger)
exp = "fixed_range"
#exp = "flexible_range"
AeolusProductObject.pickle_data_file=os.path.join("SeaFlectProjectData","DataPool","seaflect_data_"+exp + ".pkl")
figure_path = os.path.join("project_figures", "notebook", exp)
if not os.path.exists(figure_path):
os.makedirs(figure_path)
[data_pool, available_periods, available_regions] = AeolusProductObject.get_data()
AeolusProductObject.data_pool=data_pool
# now that we have a consistent dataset we can make an annual
AeolusProductObject._create_annual_(periods=available_periods)
# and globa
AeolusProductObject._create_global_(regions=available_regions)
#definitions of range bins accoring to c-convention
range_bin={"top":0, "atmosphere":3, "surface":4, "ocean":5, "background":6}
Conceptual ModelΒΆ
To understand the radiative properties of the composite layer, a simple conceptual model has been developed. This conceptual model considers the macro physical properties of the layer, and link the incoming radiation at the upper and lower interfaces of the layer to the the emerging radiation at these interfaces.
This notebook space here is too short to provide all of the equations, but in words the conceptual model indicates that the upwelling radiation at the top, is the sum of the upwelling radiation at the bottom of the layer multiplied by the transmission of the layer, and the downwelling radiation at the top of the layer multiplied by the reflection of the layer, were we ignore any internal source of radiation. As an example the upwelling radiation at the upper interface of a layer:
where the upper script gives the direction + for upwelling and - for down welling , and the sub-script the location 0 the lower interface, and 1 the upper interface.
A similar equation can be written for the emerging radiation at the lower interface.
The equations are prepared with the assumption that these macro radiative properties \(R\) and \(T\) can be linked to the micro physical properties of a homogeneous layer (e.g. Rayleigh scattering).
With the notion that the composite range bin is composed of a homogeneous atmospheric, surface and oceanograpical layer, we are interested in the macro radiative properties of the composite layer. This can be derived on the basis of the above equations, assuming that there is conservation of energy at the interfaces. The net result is that the macro physical reflection (indicated as \(R_{aso}\), were the a stands for atmosphere, s for surface and o for ocean) can be expressed in terms of these properties for the individual layers.
The second term on right hand side is the multiple reflection term. What we see from this is that the reflection of the composite layer is the reflection of the combined surface and ocean layer (which is what Li et al in their paper considers), extended with a multiple scattering term, as the radiation emerging from the surface-ocean composite layer needs to travel through the atmosphere.
In the end we are interested in an estimation of the \(R_{so}\) term, which can only be derived if we have an estimation of \(R_a\) and \(T_a\).
We can estimate \(R_a\) through the notion that in the clear atmosphere the return signal originates from Rayleigh scattering, which is a function of the number of molucules. We can estimate the Rayleigh scattering from the signal in the range bin above the surface range bin, by taking the number of molecules, into account. Which means we need to normalise signal 1 by the layer thickness in pressure space and multiply this by the pressure difference between the surface pressure and pressure level of the first interface. As we do not have the pressure values, we take the altitudes.
Thus:
were \(S\) represent the measured signal, \(\Delta z\) is the geometrical thickness, and \(z\) is the altitude, where it is acknowledged that the return signal is not proportional to geometrical thickness, but to the number of molecules in a layer.
Here we assume a cloud and aerosol free atmospheric layer. It is thus assumed that the transmission of the atmospheric layer is determined by the Rayleigh scattering only. We have no real estimation for this yet (we are building this model), we can fix it to a specific value. (0.9 in the calculations presented here)
In the end we are interested in the \(R_{so}\) as a function of the observatables. Inversion of the conceptual model indicate
it should be noted that \( \left( R_{aso} -R_{a} \right) \) is the estimation of \(R_{so}\) used by Li et al (2010). The second term (\(\left( T^2_a + R_a\left[R_{aso} - R_{a}\right] \right)\)) is a multiple scattering term which accounts for the transmission through the atmospheric layer.
Unfortunately we would need access to the incident radiation for each layer as the reflectance can be estimated as
and we measure the return signal of the layer \(F^+_1\), hence if the incident radiation is known then we can calculate the reflectance, using the auxilliary information to estimate the transmission would close the equations.
The incident radiation at the top of the surface layer is not known, and hence a proxy value for is selected. The uv_total_energy with a value of 60 is a factor 1000 too small as experiments shown. Hence a proxy value of 60000 is used in some cases, but this has no value other than bringing results within a reasonable scale.
region="E"
graphic=[]
period="AnnualCycle"
dataset=data_pool[period][region]
time = dataset["L1Bmie"]["date"]
p=figure(x_axis_type="datetime", width=800, height=150, \
tools="box_zoom,pan,crosshair, reset,save,wheel_zoom", \
title="Raso {} {}".format(period, region))
p.circle(time, dataset["Product"]["Reflectance"]["Raso"][:], color="black", \
legend_label="Composite")
p.yaxis.axis_label = 'Raso [-]'
p.y_range=Range1d(-0.05, 0.2)
p.legend.location = "top_left"
graphic.append(p)
t=figure(x_axis_type="datetime", width=800, height=150, \
tools="box_zoom,pan,crosshair, reset,save,wheel_zoom", \
title="Ra {} {}".format(period, region))
t.circle(time, dataset["Product"]["Ra"][:], color="red", legend_label="Atmosphere")
t.yaxis.axis_label = 'Ra [-]'
t.y_range=Range1d(-0.05, 0.2)
t.legend.location = "top_left"
graphic.append(t)
plot = gridplot(graphic, ncols=1, width=800, height=300, toolbar_location="right")
show(plot)
if export:
figure_name = os.path.join(figure_path,"reflectance_signal_{}_{}.png".format(period, region))
export_png(plot, filename=figure_name)
DiscussionΒΆ
The data presented above are the reflection of the composite layer (\(R_{aso}\)), and the atmospheric contribution in the surface range bin, estimated from the thickness of this atmospheric layer and the reflection of the atmospheric range bin above the surface range bin, normalised by its thickness.
Clearly seen is the effect of the geometrical thickness of the atmospheric sub-layer on the \(R_a\) signal which has large values in November 20119 and in the May - September 2020 period.
Reflection according to the simplified RT modelsΒΆ
The reflectance of the combined surface and ocean layer after removing the atmospheric contribution from the \(R_{aso}\) is shown in the next two figures. The top panel shows results from the conceptual model, while the bottom results show the results from the model adopted by Li et al. The difference between the two models is the treatment of the transmission - reflection by the atmospheric sub layer. Although the multiple scattering contribution is expected to be small, the transmission correction is significant. The transmission should be related to the number of molecules in a layer, this is currently not implemented. This could be a source of uncertainty.
region="E"
graphic=[]
period="AnnualCycle"
dataset=data_pool[period][region]
time = dataset["L1Bmie"]["date"]
p=figure(x_axis_type="datetime", width=800, height=150, \
tools="box_zoom,pan,crosshair, reset,save,wheel_zoom", \
title="Rso {} {}".format(period, region))
p.circle(time, dataset["Product"]["Reflectance"]["Cm"][:], color="green", legend_label="CM")
p.yaxis.axis_label = 'Reflectance [-]'
p.y_range=Range1d(-0.05, 0.2)
p.legend.location = "top_left"
graphic.append(p)
p=figure(x_axis_type="datetime", width=800, height=150, \
tools="box_zoom,pan,crosshair, reset,save,wheel_zoom", \
title="Rso {} {}".format(period, region))
p.circle(time, dataset["Product"]["Reflectance"]["Li"][:], color="blue", legend_label="Li")
p.yaxis.axis_label = 'Reflectance [-]'
p.y_range=Range1d(-0.05, 0.2)
p.legend.location = "top_left"
graphic.append(p)
plot = gridplot(graphic, ncols=1, width=800, height=300, toolbar_location="right")
show(plot)
if export:
figure_name = os.path.join(figure_path,"reflectance_signal_{}_{}.png".format(period, region))
export_png(plot, filename=figure_name)
Correlation with windspeedΒΆ
Finally a possible correlation of the proxy reflectance with windspeed can be presented. This is only meaningful for the observations which passed successfully the quality control. An example is shown next. This is the case for the rigorous setting of the QC (filters adopted are only positive signals, scattering ratio below threshold, SNR between expected range and standard deviation of the signal below a threshold. Values of the thresholds are presented in previous pages.) More examples are presented in the discussion section. Also shown are the theoretical results for high and low concentrations of chlorophyl.
region="E"
period="AnnualCycle"
qc_level = "sgnl_scat_snr_stdev"
#
qc = AeolusProductObject.apply_quality_control(level=qc_level, period=period, region=region)
dataset =AeolusProductObject.data_pool[period][region]
if qc == None:
time_range=slice(kwargs["time_range"][0], kwargs["time_range"][1])
else:
time_range=qc
time = dataset["L1Bmie"]["date"][time_range]
x = dataset["MET"]["surface_wind_speed"][time_range]
#basic signal
col = dataset["L1Bmie"]["snr"][time_range,AeolusProductObject.layer_index["surface"]]
# mapper = linear_cmap(field_name='col', palette=Spectral11 ,low=min(col) ,high=max(col))
mapper = linear_cmap(field_name='col', palette=Spectral11, \
low=AeolusProductObject.snr_low ,high=AeolusProductObject.snr_up)
# cm results
y = dataset["Product"]["Reflectance"]["Cm"][time_range]
source = ColumnDataSource({'x':x,'y':y,'col':col})
graphic=[]
p=figure(plot_width=400, plot_height=300, \
tools="box_zoom,pan,crosshair, reset,save,wheel_zoom", \
title="{} {} {} {}".format(period, region, qc_level, len(time)))
p.circle( "x", "y", source=source ,line_color=mapper,color=mapper,size = 20)
p.line(wind_speed, Cm["high"][0,:],color="blue", legend_label="high", line_width=8)
p.line(wind_speed, Cm["low"][0,:],color="green", legend_label="low", line_width=8)
p.xaxis.axis_label = 'wind speed [m/s]'
p.yaxis.axis_label = 'Proxy Reflectance [-]'
p.x_range = Range1d(0, 25)
p.y_range=Range1d(0,0.05)
bar = ColorBar(color_mapper=mapper['transform'], width=8, location=(0,0))
p.add_layout(bar, "left")
graphic.append(p)
plot = gridplot(graphic, ncols=1, width=800, height=500, toolbar_location="right")
show(plot)
if export:
figure_name = os.path.join(figure_path,"correlation_accd_windspeed_{}_{}.png".format(period, region))
export_png(plot, filename=figure_name)